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ABSTRACT 

Using  the  PW91  and  PBE  density  functional  theories 
(DFT),  we  have  studied  four  energetic  molecular  crystals: 
nitromethane,  HMX,  RDX,  and  CF20  with  a  wide  range 
of  basis  sets.  Our  goal  is  to  assess  the  accuracy  of  DFT 
when  applied  to  organic  molecular  crystals  (such  as 
energetic  materials)  as  scientists  are  beginning  to  include 
this  methodology  in  energetic  materials  research  without 
knowledge  of  the  limitations  of  the  method. 
Intramolecular  distances,  simple  angles,  and  band  gaps 
are  converged  at  plane  wave  cutoff  energies  (Ecut)  of  430 
to  495  eV.  Cell  parameters  were  determined  over  a  range 
of  Ecut  values  from  280  eV  to  700  or  800  eV,  depending 
upon  the  system.  Fattice  vectors,  however,  display  large 
errors  in  the  range  of  0.2  A  to  1 .0  A  (up  to  a  9.6%  error  at 
the  highest  Ecut  used  here),  and  a  very  slow  convergence 
on  basis  set  size.  We  hypothesize  the  error  in  the  lattice 
vectors  is  due  to  a  lack  of  van  der  Waals  forces  in  current 
DFT  functionals.  This  deficiency  will  have  unforeseen 
consequences  on  all  crystal  calculations  for  organic 
molecules,  and  therefore  caution  should  be  employed 
whenever  interpreting  results  obtained  from  the  current 
DFT  functionals  available  in  solid  state  codes.  To 
properly  describe  the  electronic  structure  of  these  types  of 
crystals,  these  results  suggest  the  need  for  new  methods 
involving  DFT  to  be  developed  which  include  accurate 
dispersion  forces. 

1.  INTRODUCTION 

Recent  advances  in  density  functional  theory  (DFT) 
(Hohenberg  and  Kohn,  1964)  have  allowed  the  study  of 
large  polyatomic  molecules  and  crystals.  These 
advancements  have  become  part  of  the  standard  set  of 
tools  used  in  the  development  and  evaluation  of  energetic 
materials.  Density  functional  theories  exhibit  reasonable 
accuracies  coupled  with  computational  efficiency,  and  for 
such  reasons  have  been  used  to  study,  both  statically  and 
dynamically,  large  molecular  systems.  Previous  work  for 
energetic  species  utilizing  density  functional  theories  has 
ranged  from  gas  phase  structure  studies  of  RDX  (Rice  and 
Chabalowski,  1997)  to  shocked  nitromethane  calculations 
(Reed  et  al.,  2000;  Smith  and  Bharadwaj,  1999).  Density 
functional  theory  has  also  been  used  in  the  charting  of  the 
unimolecular  decomposition  pathways  for  both  RDX 


(Chakraborty  et  al.,  2000)  and  HMX  (Chakraborty  et  al., 

2001). 

Inspired  by  the  work  of  Reed  et  al.  (Reed  et  al., 
2000),  we  initially  desired  to  examine  the  possibility  of 
bandgap  lowering  under  compression  as  a  mechanism  for 
detonation,  as  proposed  by  Kuklja  et  al  (Kuklja  and  Kunz, 
1999).  However,  upon  testing  for  basis  set  convergence, 
we  realized  there  were  significant  errors  present  in  the 
lattice  vectors  when  using  density  functional  theories. 
These  errors  in  the  lattice  vectors  can  translate  into  large 
errors  in  the  computed  crystalline  densities  of  these 
energetic  molecular  crystals. 

The  density  of  crystalline  energetic  materials  is 
considered  to  be  one  of  the  primary  physical  parameters 
in  detonation  performance,  since  detonation  velocities  and 
pressures  of  energetic  materials  are  proportional  to 
powers  in  the  density.  Currently,  most  density  predictions 
are  within  2-4%  of  experimental  values,  depending  on  the 
class  of  compounds.  Experimentalists  would  prefer 
predictions  with  deviations  of  1%  or  less  (Dr.  Rao 
Surapaneni,  ARDEC,  Picatinny  Arsenal,  NJ,  private 
communication).  Thus,  any  calculations  that  yield 
computed  densities  with  errors  greater  than  4-5%  are 
undesirable.  In  this  paper,  we  shall  demonstrate  that 
current  formulations  of  DFT  are  unsuitable  for  predictions 
of  densities  of  energetic  molecular  crystals  at  ambient 
conditions,  and  shall  quantify  the  errors  one  can  expect 
from  density  functional  theory. 

2.  COMPUTATIONAL  DETAILS 

Calculations  of  the  crystal  parameters  and 
intramolecular  bond  distances  and  angles  for 
nitromethane  (N02CH3),  1 ,3,5,7-tetranitro- 1 ,3,5,7-tetra- 
azacyclo-octane  (HMX),  cyclotrimethylene-trinitramine 
(RDX),  and  2,4,6,8,10,12-hexanitrohexa-azaisowurzitane 
(CL20)  were  done  using  the  generalized  gradient 
approximation  density  functional  theory  (GGA  DFT) 
Perdew-Wang  91  (PW91)  (Perdew  et  al.,  1991)  functional 
as  implemented  in  the  Vienna  Ab-Initio  Package  (VASP) 
(Kresse  and  Furthmuller,  2003).  The  Perdew-Burke- 
Emzerhof  (PBE)  (Perdew  et  al.,  1996)  functional  was  also 
employed  as  an  alternate  method  for  nitromethane  and 
RDX.  The  Vanderbilt  ultrasoft  pseudopotentials  (USP) 
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(Vanderbilt,  1990)  and  Monkhorst-Pack  k-point 
generation  method  were  used  for  all  calculations.  Plane 
wave  kinetic  energy  cut-offs  (Ecut)  ranged  from  280  eV  to 
800  eV. 

Starting  from  experimental  crystal  structures,  all 
ionic  and  cell  degrees  of  freedom  were  allowed  to  relax 
concurrently  and  without  enforcing  symmetry  restrictions. 
Theoretical  equilibrium  crystal  structures  were  obtained 
through  analytic  gradient  techniques.  The  electronic 
energies  were  converged  to  2.0x1 0"6  eV,  while  the 
structures  were  converged  once  the  difference  in  free 
energy  between  gradient  steps  was  less  than  2.0x1 0"5  eV. 
As  these  molecular  crystals  are  insulators,  setting  the 
smearing  width  to  0.0001  in  the  Methfessel-Paxton 
method  (Methfessel  and  Paxton,  1989)  minimized 
electron  smearing.  To  accelerate  convergence,  the 
RMM-DIIS  method  was  employed,  except  for  certain 
calculations  where  it  was  necessary  to  use  conjugate 
gradients  instead.  As  a  comparison  of  plane  wave  versus 
Gaussian  basis  set  calculations,  we  obtained  calculations 
of  nitromethane  using  the  Guassian03  (G03)  solid-state 
package.  (Frisch  et  al.,  2003)  Two  calculations  were 
performed  using  the  3-21G  and  6-3 1G**  basis  sets. 
(Binkley  et  ah,  1980;  Gordon  et  ah,  1983;  Hariharan  and 
Pople,  1973;  Francl  et  ah,  1982) 

3.  RESULTS  AND  DISCUSSION 
3.1  Nitromethane 

The  simplest  of  the  four  energetic  materials  studied  is 
nitromethane,  with  four  molecules  contained  in  the 
primitive  cell.  The  space  group  is  P2!2i2i  with  28  atoms 
in  the  cell  and  experimental  lattice  vectors  (taken  at  4.2 
K)  of  a  =  5.183  A,  b  =  6.236  A,  and  c  =  8.518  A  and 
lattice  angles  of  a  =  P  =  y  =  90.0°.  (Trevino  et  ah,  1980) 
We  tested  for  an  appropriate  k-point  grid  by  running 
calculations  at  a  kinetic  energy  cut-off  of  545  eV  with  one 
(lxlxl),  four  (2x2x2),  fourteen  (3x3x3),  and  thirty-two 
(4x4x4)  k-points  with  the  PW91  functional.  The 
difference  in  lattice  vectors  between  the  four  and  thirty- 
two  k-point  grids  was  at  most  -0.007  A  (-0.118%). 
Therefore,  the  remainder  of  the  calculations  were 
performed  using  a  four  k-points  grid.  In  testing  for 
convergence  with  kinetic  energy  cut-off,  we  ran 
calculations  at  280,  330,  380,  430,  495,  545,  645,  and  700 
eV  with  the  PW91  functional. 

Inspecting  first  the  intramolecular  bond  distance 
errors  versus  experiment  in  Table  1,  we  see  by  an  Ecut  of 
495  eV,  the  errors  have  converged  to  very  acceptable 
values.  With  an  absolute  mean  of  only  0.017  A  and  a 
standard  deviation  of  0.013  A  by  495  eV,  we  see  an 
excellent  agreement  with  experiment.  The  largest  error 
for  nitromethane  is  one  of  the  nitrogen-oxygen  bonds, 
with  every  kinetic  energy  cut-off  consistently  resulting  in 


a  0.04  A  error.  The  PW91  functional  yields 
approximately  the  same  bond  lengths  for  each  of 
nitrogen-oxygen  bonds,  while  the  experimental  lengths  of 
the  N-0  bond  differ  (1.209  A  and  1.223  A).  This  error  in 
the  N-0  bond  length  is  not  carried  over  into  the  angles 
(Table  2),  where  it  is  one  of  the  methyl  hydrogens  that 
yield  the  largest  error  of  approximately  3.6  degrees 
relative  to  the  nitrogen.  Outside  of  this  lone  case,  no 
other  angles  exhibits  any  error  to  experiment  greater  than 
3  degrees  and  are  for  the  most  part  well  behaved. 

Table  1:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  bond 
distances  versus  experiment  for  various  kinetic  energy 
cut-offs  for  the  nitromethane  crystal  with  PW91.  _ 


Ecut 

(eV) 

Mean 

(A) 

% 

Mean 

Abs 

Mean 

(A) 

Std 

Dev 

(A) 

Max 

Error 

(A) 

280 

0.0063 

0.0630 

0.0190 

0.0231 

0.037 

330 

0.0115 

1.0309 

0.0172 

0.0182 

0.038 

380 

0.0135 

1.1824 

0.0158 

0.0155 

0.037 

430 

0.0155 

1.3325 

0.0158 

0.0136 

0.037 

495 

0.0167 

1.4239 

0.0167 

0.0128 

0.037 

545 

0.0167 

1.4239 

0.0167 

0.0128 

0.037 

645 

0.0168 

1.4352 

0.0168 

0.0126 

0.037 

700 

0.0167 

1.4237 

0.0167 

0.0125 

0.037 

Table  2:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  angles 
versus  experiment  for  various  kinetic  energy  cut-offs  for 
the  nitromethane  crystal  with  PW9 1 .  _ _ 


Ecut 

(eV) 

Mean 

(°) 

% 

Mean 

Abs 

Mean 

c°) 

Std 

Dev  (°) 

Max 

Error 

n 

280 

0.0333 

0.0548 

1.5444 

2.1319 

4.5 

330 

0.0222 

0.0422 

1.3556 

1.9942 

4.2 

380 

0.0111 

0.0318 

1.3000 

1.8878 

3.9 

430 

0.0111 

0.0299 

1.3000 

1.8128 

3.7 

495 

0.0111 

0.0300 

1.3000 

1.8114 

3.7 

545 

0.0222 

0.0376 

1.3111 

1.7513 

3.4 

645 

0.0333 

0.0484 

1.3444 

1.8104 

3.6 

700 

0.0222 

0.0392 

1.1068 

1.7957 

3.6 

Turning  now  to  the  crystal  lattice  parameters,  we  see 
a  disturbing  trend  in  the  error  to  experiment  versus  the 
kinetic  energy  cut-off.  Figure  1  illustrates  how  at  high 
Ecut,  the  percent  error  to  experiment  is  fairly  large,  with  a 
maximum  percent  error  of  8.  5  %  (0.53  A)  for  the  b  lattice 
vector  at  545  eV.  While  the  negative  errors  seen  for  the 
small  Ecut  are  also  large,  this  is  understandable  as  there  are 
far  too  few  plane  waves  at  these  low  Ecut  to  properly 
describe  the  momentum  space.  However,  at  larger  Ecut 
values,  the  convergence  of  the  b  lattice  vector  to 
approximately  7%  error  even  at  a  large  Ecut  is  indicative 
that  the  current  implementation  of  DFT  is  unsuited  to 
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molecular  crystalline  compounds.  We  shall  see  similar 
trends  in  the  other  molecular  crystals  in  this  study. 

A  possible  source  of  error  for  these  calculations  is 
that  VASP  keeps  the  number  of  plane  waves  in  the 
calculation  constant,  regardless  of  change  in  cell  volume. 
To  be  certain  that  we  properly  saturated  the  momentum 
space  with  plane  waves  and  that  increasing  the  cell 
volume  will  not  cause  a  Pulay  stress  error,  we  started  a 
new  cell  convergence  calculation  from  the  converged 
geometry  at  545  eV  without  changing  any  parameters. 
Upon  completion,  the  largest  change  in  lattice  vector  was 
-0.024  A  (-0.47  %)  for  the  a  lattice  vector.  This  minimal 
change  allows  us  to  state  that  at  the  higher  Ecut  the 
momentum  space  is  adequately  described  and  that  Pulay 
stresses  are  not  a  concern. 


Figure  1:  PW91  Percent  Errors  versus 
Experimental  Lattice  Vectors  for  Nitromethane 


Kinetic  Energy  Cut-off  (eV) 


Figure  2:  Percent  Errors  versus  Experimental 
Lattice  Vectors  for  Nitromethane 


Lattice  Vector 

Other  possible  sources  of  error  would  be  errors 
inherent  to  the  PW91  functional  or  the  plane  wave  basis 
sets.  We  therefore  computed  the  crystalline  dimensions 
using  the  VASP  implemented  PBE  functional  with  a  four 
k-point  grid  and  with  the  495,  545,  and  645  eV  basis  sets. 
We  also  determined  the  lattice  vectors  using  the  3-2 1G 


and  6-3 1G**  Gaussian  basis  sets  with  the  G03 
implemented  PBE  functional.  Figure  2  illustrates  the 
errors  in  the  a,  b,  and  c  lattice  vectors  for  the  above 
methods  and  functionals,  as  well  as  with  the  PW91 
functional  using  the  495,  545,  and  645  eV  basis  sets  for 
comparison. 

The  plane  wave  PBE  functional  yields  errors 
approximately  0.5- 1.0%  lower  than  the  PW91  functional, 
with  a  few  exceptions.  For  the  b  lattice  vector,  PW91  at 
545  eV  has  larger  errors  than  the  PBE  functional,  while 
the  reverse  holds  true  for  the  645  eV  basis  set.  The  PBE 
functional  does  not  seem  to  dramatically  improve  upon 
the  results  found  with  the  PW91  functional. 

For  the  Gaussian  basis  sets,  it  is  readily  apparent  the 
3-21G  basis  set  is  too  small  to  properly  describe  the 
nitromethane  crystal.  The  6-3 1G**  basis  set,  however, 
does  reasonably  well  for  the  a  and  c  lattice  vectors,  but 
shows  a  4.96%  error  in  the  b  lattice  length.  It  seems  the 
Gaussian  basis  sets  display  similar  trends  in  the  errors  as 
the  plane  wave  basis  sets,  and  as  the  results  should  be 
independent  of  the  nature  of  the  basis  set  employed,  this 
was  not  unexpected. 

3.2  HMX 

We  look  next  to  P-HMX,  the  most  stable  polymorph 
of  HMX,  which  exhibits  a  space  group  of  P  2x!z  with  56 
atoms  in  the  unit  cell  and  experimental  lattice  vectors  of  a 
=  6.54  A,  b  =  11.05  A,  and  c  =  8.7  A  and  a  lattice  angles 
of  a  =  y  =  90.0°,  and  p  =  124.3°.  (Choi  and  Boutin,  1970) 
We  tested  for  an  appropriate  k-point  grid  by  running 
calculations  at  a  kinetic  energy  cut-off  of  545  eV  with  two 
(2x1x2)  and  thirty- two  (4x4x4)  k-points.  The  largest 
difference  in  lattice  vectors  between  the  two  different  k- 
point  grids  was  0.009  A  (0.078%).  Therefore  the 
calculations  for  this  study  were  performed  using  the  two 
(2x1x2)  k-points  grid.  In  testing  for  convergence  with 
kinetic  energy  cut-off,  we  ran  calculations  at  280,  330, 
380,  430,  495,  545,  and  700  eV. 

As  one  can  see  in  Tables  3  and  4,  HMX  follows  a 
similar  trend  as  nitromethane,  where  by  an  Ecut  of  495  eV 
the  errors  have  converged.  We  see  a  convergence  to  an 
absolute  mean  error  of  0.014  A  (0.8  %  mean  error)  in 
bond  lengths  and  a  standard  deviation  convergence  to 
0.015  A  compared  to  experiment.  The  bond  that  displays 
the  largest  theoretical  error  is  a  nitrogen-oxygen  bond  off 
an  axial  N02  group.  However,  both  N-0  bonds  of  the 
axial  N02  groups  and  one  N-0  bond  of  the  equatorial 
N02  groups  display  errors  of  approximately  0.03  A.  No 
other  bond  exhibits  errors  larger  than  0.019  A  in  any  Ecut 
after  280  eV,  which,  as  stated  in  the  previous  section,  is 
an  unsuitable  Ecut  for  molecular  crystal  calculations.  All 
bond  angles  are  within  ±  3°  of  experiment  for  all  values 
of  Ecut  (Table  4). 
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Table  3:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  bond 
distances  versus  experiment  for  various  kinetic  energy 
cut-offs  for  the  P-HMX  crystal  with  PW91. 


Ecut 

(eV) 

Mean 

(A) 

% 

Mean 

Abs 

Mean 

(A) 

Std 

Dev 

(A) 

Max 

Error 

(A) 

280 

-0.0031 

-0.1911 

0.0195 

0.0229 

0.045 

330 

0.0035 

0.2973 

0.0135 

0.0190 

0.044 

380 

0.0074 

0.5782 

0.0116 

0.0164 

0.041 

430 

0.0091 

0.7080 

0.0125 

0.0156 

0.039 

495 

0.0104 

0.7994 

0.0136 

0.0154 

0.038 

545 

0.0104 

0.8056 

0.0137 

0.0153 

0.037 

700 

0.0104 

0.8047 

0.0137 

0.0154 

0.038 

Table  4:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  angles 
versus  experiment  for  various  kinetic  energy  cut-offs  for 
the  P-HMX  crystal  with  PW9 1 .  _ _ 


Ecut 

(eV) 

Mean 

(°) 

% 

Mean 

Abs 

Mean 

c°) 

Std 

Dev  (°) 

Max 

Error 

o 

280 

-0.0400 

-0.0233 

0.8880 

1.3150 

3.3 

330 

0.0760 

0.0763 

0.7400 

1.1107 

2.6 

380 

0.0640 

0.0663 

0.7280 

1.0618 

2.6 

430 

0.0600 

0.0600 

0.7960 

1.0851 

-2.6 

495 

0.0760 

0.0731 

0.8520 

1.1092 

-2.7 

545 

0.0800 

0.0764 

0.8720 

1.1273 

-2.8 

700 

0.0520 

0.0523 

0.8200 

1.0978 

-2.8 

Figure  3:  PW91  Percent  Errors  versus 
Experimental  Lattice  Vectors  for  HMX 


280  330  380  430  495  545  700 
Kinetic  Energy  Cut-off  (eV) 


In  Figure  3,  one  sees  that  the  percent  errors  in  the 
lattice  vectors  increases  as  the  Ecut  values  increase. 
Similarly  to  nitromethane,  below  430  eV  the  lattice 
vectors  are  all  shorter  than  experiment,  while  for  430  eV 
and  above  the  percent  errors  were  positive  and  growing 
with  increasing  Ecut.  However,  unlike  nitromethane,  all 
the  lattice  vectors  appear  to  begin  converging  by  495  eV. 
The  percent  errors  are  also  not  as  large  as  those  seen  in 
nitromethane,  ranging  from  approximately  3%  for  the  c 


vector  to  5%  for  the  a  vector  at  an  Ecut  of  700  eV.  The 
percent  error  in  the  p  lattice  angle  shows  no  pattern  and  is 
consistently  small,  with  most  errors  under  0.34°,  with  the 
exception  of  -1.5°  at  545  eV.  For  this  molecular  crystal, 
both  the  intramolecular  errors  and  the  lattice  vector  errors 
seem  to  be  converged  by  495  eV.  We  shall  see  that  these 
are  the  best  results  for  DFT  within  the  sample  set  of  four 
molecular  crystals  in  this  study. 

3.3  RDX 

Examining  next  RDX,  we  are  fortunate  enough  have 
experimental  gas  phase  data  to  compare  against.  In 
comparing  to  gas  phase  data,  we  will  be  able  to  gauge  the 
accuracy  of  the  theory  with  respect  to  plane  waves  and 
pseudopotentials.  As  there  also  exist  previous  theoretical 
calculations  on  the  gas  phase  data  using  atomic  orbitals, 
we  will  compare  to  these  results  as  well. 

As  VASP  was  designed  to  compute  periodic  systems, 
to  approximate  gas  phase  conditions  we  ran  an  “isolated” 
molecule  in  a  15x15x15  A3  box  with  one  k-point. 
Although  this  cell  is  replicated  in  all  three  spatial 
dimensions,  the  large  vacuum  around  each  molecule 
should  be  sufficient  to  eliminate  any  forces  felt  from  its 
neighbors.  Calculations  were  done  using  kinetic  energy 
cut-offs  of  280,  380,  and  495  eV  and  all  calculations  were 
started  using  the  gas  phase  experimental  structure  as  the 
initial  geometry. 

Table  5:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  bond  distances  versus 
experiment  for  PW91,  MP2,  and  B3LYP  for  the  gas  phase 
RDX  molecule. 


Mean 

(A) 

% 

Mean 

Abs 

Mean 

(A) 

Std 

Dev 

(A) 

Max 

Error 

(A) 

PW91 
(280  eV) 

0.0279 

2.2232 

0.0282 

0.0220 

0.070 

PW91 
(380  eV) 

0.0186 

1.5197 

0.0205 

0.0171 

0.045 

PW91 
(495  eV) 

0.0181 

1.4860 

0.0201 

0.0171 

0.045 

MP2/ 

6-3 1G* 

0.0184 

1.4673 

0.0197 

0.0180 

0.060 

B3LYP/ 

6-3 1G* 

0.0146 

1.1558 

0.0166 

0.0160 

0.051 

B3LYP/6- 

31 1+G** 

0.0128 

0.9965 

0.0156 

0.0171 

0.054 

Comparing  to  the  theoretical  work  done  previously 
by  Rice  and  Chabalowski  (Tables  5  and  6)  (Rice  and 
Chabalowski,  1997),  we  see  that  by  380  eV,  PW91  with 
plane  waves  yields  similar  results  to  MP2/6-31G*.  There 
is  little  improvement  in  going  to  a  larger  Ecut,  as  the 
results  have  converged  by  380  eV.  B3LYP  with  a 
Gaussian  basis  set,  however,  displays  a  0.0035  A  decrease 
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in  absolute  mean  error  with  the  6-3 1G*  basis  set  versus 
PW91  with  a  495  eV  Ecut.  There  is  marginal 
improvement  in  the  errors  for  B3LYP  in  going  to  the 
larger  6-311+G**  basis  set,  decreasing  the  absolute  mean 
error  by  only  an  additional  0.001  A.  The  bond  with  the 
largest  error  for  all  methods  and  basis  sets  compared  to 
gas  phase  data  (except  for  PW91  at  280  eV)  is  the 
equatorial  nitrogen-nitrogen  bond.  The  accuracy  of  these 
results  validates  the  use  of  ultrasoft  pseudopotentials,  as 
they  are  sufficient  to  describe  the  intramolecular  structure 
of  organic  species.  These  results  also  indicate  that  for 
accurate  structures  one  does  not  need  very  high  kinetic 
energy  cut-offs  for  intramolecular  distances.  We  have 
already  seen  evidence  of  this  for  the  nitromethane  and 
HMX  molecular  crystals.  As  we  shall  continue  to 
demonstrate,  however,  the  crystal  lattice  vectors  are  very 
sensitive  to  the  level  of  Ecut  chosen. 

Table  6:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  angles 
versus  experiment  for  PW91,  MP2,  and  B3LYP  for  the 
gas  phase  RDX  molecule. _ _ _ _ 


Mean  (°) 

% 

Mean 

Abs 

Mean 

(°> 

Std 

Dev 

(°) 

Max 

Error 

h 

PW91 
(280  eV) 

-0.0306 

-0.0306 

0.9750 

1.1181 

2.2 

PW91 
(380  eV) 

0.2056 

0.1739 

0.8722 

1.0450 

2.3 

PW91 
(495  eV) 

0.2194 

0.1863 

0.9306 

1.0828 

2.2 

MP2/ 

6-3 1G* 

-0.8367 

-0.7110 

0.1459 

1.9304 

-7.25 

B3LYP/ 

6-3 1G* 

-0.2794 

-0.2352 

1.0506 

1.5470 

-5.81 

B3LYP/6- 

311+G** 

0.2292 

-0.1928 

0.9897 

1.4470 

-5.48 

ran  calculations  at  280,  330,  380,  430,  495,  545,  700,  and 
800  eV. 

In  general  the  statistical  errors  versus  experiment  for 
RDX  are  slightly  better  in  the  crystal  phase  than  in  the  gas 
phase.  Absolute  means  have  converged  to  0.018  A  (1.4 
mean  percent  error)  in  the  crystal  phase,  an  improvement 
of  0.002  A  over  the  gas  phase  predictions  (Table  7).  The 
bond  with  the  largest  error  varies  with  Ecut  for  this 
molecular  crystal.  For  330  and  380  eV,  an  axial  N-0  is 
the  longest  bond,  while  for  430  through  800  eV  the  short 
C-H  bond  (1.059  A  in  experiment)  yields  the  largest  error 
versus  experiment,  as  theory  predicts  all  C-H  bonds  to  be 
approximately  the  same  length  (-1.09  A).  Generally,  all 
N-0  bonds  are  too  long  and  have  errors  on  the  order  of 
0.04  A,  whether  they  are  axial  or  equatorial  N02  groups. 
For  Ecut  larger  than  380  eV,  almost  all  bonds  are  longer 
than  experimental  values  and  have  converged  to  the  same 
value  by  430  eV. 

Table  7:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  bond 
distances  versus  experiment  for  various  kinetic  energy 
cut-offs  for  the  RDX  crystal  with  PW91. _ _ 


Ecut 

(eV) 

Mean 

(A) 

% 

Mean 

Abs 

Mean 

(A) 

Std 

Dev 

(A) 

Max 

Error 

(A) 

280 

-0.0058 

-0.3029 

0.0256 

0.0292 

-0.047 

330 

0.0058 

0.5679 

0.0173 

0.0217 

0.041 

380 

0.0139 

1.1729 

0.0155 

0.0160 

0.041 

430 

0.0158 

1.309 

0.0170 

0.0146 

0.040 

495 

0.0166 

1.3702 

0.0175 

0.0144 

0.040 

545 

0.0167 

1.3780 

0.0178 

0.0145 

0.041 

700 

0.0165 

1.3619 

0.0175 

0.0148 

0.040 

800 

0.0167 

1.3794 

0.0177 

0.0146 

0.041 

We  examine  next  the  crystal  phase  of  a-RDX,  the 
room  temperature  polymorph  of  crystalline  RDX, 
containing  168  atoms  per  cell  and  with  a  space  group  of  P 
bca  and  experimental  lattice  vectors  of  a  =  13.182  A,  b  = 
1 1.574  A,  and  c  =  10.709  A  and  lattice  angles  of  a  =  P  =  y 
=  90.0°  (Choi  and  Prince,  1972).  K-point  testing  was 
performed  using  the  one  (lxlxl)  and  four  (2x2x2)  k-point 
grids  at  495  eV.  The  largest  absolute  difference  in  lattice 
vectors  between  the  two  different  k-point  grids  was  the  c 
vector  at  0.15  A  (1.4%)  and  the  smallest  difference  the  a 
vector  at  0.02  A  (0.17%),  with  the  four  k-point  grid 
yielding  a  longer  a  lattice  vector  but  a  shorter  c  vector 
compared  to  the  one  k-point  grid.  As  the  c  vector  still 
exhibits  a  large  percent  error  of  7.78%  with  the  four  k- 
point  grid,  coupled  with  the  size  of  unit  cell,  extensive 
calculations  with  the  four  k-point  grid  were  not 
computationally  feasible,  and  therefore  the  calculations 
were  performed  using  the  one  (lxlxl)  k-point  grid.  In 
testing  for  convergence  with  kinetic  energy  cut-off,  we 


Table  8:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  angles 
versus  experiment  for  various  kinetic  energy  cut-offs  for 
the  RDX  crystal  with  PW91. _ _ _ 


Ecut 

(eV) 

Mean 

(°> 

% 

Mean 

Abs 

Mean 

(°) 

Std 

Dev  (°) 

Max 

Error 

(°) 

280 

0.0917 

0.0898 

1.4472 

2.0136 

5.0 

330 

0.1778 

0.1582 

0.9278 

1.1931 

3.4 

380 

0.1917 

0.1662 

0.5361 

0.6487 

2.0 

430 

0.2389 

0.2066 

0.5444 

0.6750 

1.8 

495 

0.3000 

0.2595 

0.6333 

0.7878 

2.2 

545 

0.2944 

0.2542 

0.6667 

0.7950 

2.1 

700 

0.3139 

0.2711 

0.6972 

0.8153 

2.3 

800 

0.3194 

0.2764 

0.6972 

0.8605 

1.9 

Similarly,  the  bond  angles  have  roughly  converged 
by  495  eV  (Table  8).  The  angles  with  the  largest  error 
varies  as  the  basis  set  increases:  a  CNO  angle  largest  at 
330  eV,  a  HCH  angle  at  430  eV,  and  for  494  through  800 
eV  the  largest  error  is  in  a  CNN  angle.  However,  as  we 
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have  seen  before,  the  errors  in  the  angles  are  on  average 
smaller  than  0.7  degrees. 


there  still  exists  a  quite  large  approximate  8%  error  (0.8 
A),  inconsistent  with  accurate  calculations. 


Figure  4:  PW91  Percent  Errors  versus 
Experimental  Lattice  Vectors  for  RDX 
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Figure  5:  Percent  Errors  versus 
Experimental  Lattice  Vectors  for  RDX 
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Figure  4  illustrates  the  errors  in  the  crystal  lattice 
vectors  for  RDX  versus  Ecut.  As  seen  previously,  the  low 
kinetic  energy  cut-offs  are  inadequate  to  properly  describe 
the  momentum  space,  and  therefore  yield  extremely  large 
errors.  However,  unlike  HMX,  even  though  the 
intramolecular  bonds  and  angles  have  converged  to 
approximately  the  same  accuracy  by  430  eV,  the  lattice 
vectors  are  still  increasing  in  length  for  the  higher  Ecut.  At 
800  eV,  RDX  exhibits  the  largest  error  in  this  study,  with 
a  9.62  %  error  (1.03  A)  in  the  c  lattice  vector,  and  large 
errors  in  the  two  remaining  lattice  vectors  with  a  0.57  A 
and  0.24  A  errors  in  the  a  and  b  lattice  vectors, 
respectively.  These  errors  in  the  lattice  vectors  are 
unfortunately  not  limited  only  to  the  PW91  functional,  but 
are  also  present  in  the  PBE  functional.  In  Figure  5,  one 
can  observe  how  closely  the  PBE  functional  tracks  the 
same  errors  as  the  PW91  functional  at  the  430  and  495  eV 
basis  sets.  While  the  PBE  functional  does  improve 
marginally  on  the  errors  in  the  c  lattice  vector,  at  495  eV 


3.4  CL20 

Finally  turning  to  S-CL20,  the  most  stable  polymorph 
of  CL20  with  144  atoms  per  crystal  cell,  which  has  a 
space  group  of  P  2\ln  and  experimental  lattice  vectors  of  a 
=  8.852  A,  b  =  12.556  A,  and  c  =  13.386  A  and  lattice 
angles  of  a  =  y  =  90.0°,  and  (3  =  106.82°.  (Nielson  et  al., 

1998)  K-point  testing  was  done  with  the  one  (lxlxl)  and 
four  (2x2x2)  k-point  grids  at  495  eV.  The  largest 
difference  in  lattice  vectors  between  the  two  different  k- 
point  grids  was  the  b  vector  at  0.19  A  (1.51%)  and  the 
smallest  difference  the  a  vector  at  0.10  A  (1.12%),  with 
the  four  k-point  grid  yielding  the  longer  lattice  vectors. 
However,  with  such  a  large  unit  cell,  extensive 
calculations  with  the  four  k-point  grid  were 
computationally  intractable,  and  therefore  the  calculations 
were  performed  using  the  one  (lxlxl)  k-points  grid.  In 
testing  for  convergence  with  kinetic  energy  cut-off,  we 
ran  calculations  at  280,  330,  380,  430,  495,  545,  and  700 
eV. 

Table  9:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  bond 
distances  versus  experiment  for  various  kinetic  energy 
cut-offs  for  the  S-CL20  crystal  with  PW91. _ _ 


Ecut 

(eV) 

Mean 

(A) 

% 

Mean 

Abs 

Mean 

(A) 

Std 

Dev 

(A) 

Max 

Error 

(A) 

280 

-0.0288 

-2.0573 

0.0293 

0.0245 

-0.090 

330 

-0.0137 

-0.9801 

0.0148 

0.0149 

-0.047 

380 

-0.0016 

-0.1206 

0.0033 

0.0045 

-0.016 

430 

0.0013 

0.0856 

0.0026 

0.0040 

0.013 

495 

0.0024 

0.1582 

0.0027 

0.0047 

0.017 

545 

0.0022 

0.1450 

0.0029 

0.0050 

0.018 

700 

0.0025 

0.1672 

0.0030 

0.0050 

0.017 

Table  10:  Mean,  percent  mean,  absolute  mean,  standard 
deviation,  and  maximum  error  in  intramolecular  angles 
versus  experiment  for  various  kinetic  energy  cut-offs  for 
the  S-CL20  crystal  with  PW91. 


Ecut 

(eV) 

Mean 

(°) 

% 

Mean 

Abs 

Mean 

c°) 

Std 

Dev  (°) 

Max 

Error 

n 

280 

-0.0972 

-0.0648 

1.7722 

2.3761 

-8.1 

330 

0.0306 

0.0377 

1.1389 

1.4587 

5.6 

380 

-0.0347 

-0.0251 

0.4458 

0.6070 

2.0 

430 

-0.0333 

-0.0262 

0.1528 

0.2195 

-0.8 

495 

-0.0014 

-0.0011 

0.0014 

0.0118 

-0.1 

545 

0.0286 

0.2311 

0.1369 

0.2161 

1.0 

700 

-0.0069 

-0.0058 

0.0264 

0.0513 

-0.1 

For  the  CL20  crystal,  intramolecular  errors  are  small 
(Table  9),  almost  an  order  of  magnitude  smaller  than  the 


6 


other  crystals  in  this  study.  Maximum  errors  are  roughly 
a  factor  of  three  smaller  than  those  seen  for  the  other 
crystals.  One  possible  explanation  of  this  phenomenon  is 
the  more  rigid  structure  of  the  CL20  molecule.  As  all 
calculations  were  initiated  from  the  experimental  crystal 
structure,  the  lack  of  atomic  mobility  restricted  the 
molecular  crystal  to  remain  close  to  the  starting  structure. 
As  before,  intramolecular  distances  are  converged  by  495 
eV. 

Table  10  shows  that  the  bond  angles  are  also 
predicted  more  accurately  for  C120  than  those  of  the  other 
crystals  in  this  study.  We  observe  angle  errors 
approximately  two  to  ten  times  smaller  than  in  previous 
structures.  While  the  errors  are  in  general  smaller  than 
for  other  structures,  there  appears  to  be  little  convergence 
in  the  angle  errors  for  increasing  basis  set  size.  However, 
with  maximum  errors  on  the  order  of  1  degree  or  smaller, 
and  standard  deviations  of  0.22  degrees  or  less  for  all 
calculations  with  basis  sets  larger  than  430  eV,  a  small 
amount  of  fluctuation  is  not  a  serious  concern. 


of  basis  sets.  Intramolecular  distances,  angles,  and  band 
gaps  are  converged  by  Ecut  of  430  to  495  eV.  Lattice 
vectors,  however,  display  large  errors  in  the  range  of  0.2 
A  to  1.0  A  (up  to  a  9.6%  error)  and  a  slower  convergence 
on  basis  set  size.  Table  1 1  illustrates  the  impact  in  these 
large  lattice  vector  errors  on  the  density,  where  we  see 
percent  differences  versus  experiment  of  10%  or  more  in 
the  computed  crystalline  densities.  Recalling 
experimentalists  would  prefer  predictions  with  deviations 
of  1%  or  less,  these  errors  in  the  crystal  densities  establish 
that  current  density  functional  theories  are  unsuitable  for 
accurate  predictions. 


Table  11:  Percent  errors  in  the  computed  crystal  densities 
versus  experiment  for  the  nitromethane,  HMX,  RDX,  and 
CL20  crystals  with  PW91. 


Ecu, 

(eV) 

Nitromethane 

HMX 

RDX 

CL20 

495 

-13.5 

-10.6 

-13.3 

-11.8 

545 

-16.3 

-12.5 

-13.3 

-13.4 

700 

-15.1 

-10.1 

-11.7 

-11.8 

However,  for  lattice  vectors  the  trend  for  large  errors 
continues,  where  we  observe  maximum  b  lattice  vector 
error  of  6.1%  (0.77  A)  at  545  eV  (Figure  6).  Even 
progressing  to  a  700  eV  Ecut,  the  b  vector  error  remains 
large  at  5.3%  (0.66  A).  While  the  percent  errors  are  on 
par  with  those  of  nitromethane  and  HMX,  this  is 
deceiving  as  the  experimental  lattice  vectors  for  CL20  are 
larger  than  those  of  the  other  two  molecular  crystals. 
Therefore  the  absolute  errors  in  angstroms  are  larger  for 
CL20  than  for  nitromethane  and  HMX.  The  disparity  in 
the  accuracy  of  the  intramolecular  distances  and  angles 
compared  to  the  inaccuracy  of  the  lattice  vectors  is  most 
pronounced  in  CL20. 

Figure  6:  PW91  Percent  Errors  versus 
Experimental  Lattice  Vectors  for  CL20 
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4.  CONCLUSION 

Using  the  PW91  and  PBE  DFT  functionals,  we  have 
studied  four  common  energetic  molecular  crystals: 
nitromethane,  HMX,  RDX,  and  CL20  with  a  wide  range 


We  tested  different  k-point  spaces  to  determine  the 
proper  saturation  level.  We  also  compared  the  PW91  and 
PBE  functional  for  both  nitromethane  and  RDX.  There 
were  minimal  differences  between  the  two  methods  for 
the  same  basis  set  size.  Gaussian  based  DFT  calculations 
using  the  PBE  functional  were  computed  for  nitromethane 
in  order  to  test  inaccuracies  in  plane  wave  basis  sets.  The 
6-3 1G**  basis  set  yielded  results  for  cell  vectors  in 
slightly  better  agreement  with  experiment  than  the  plane 
wave  results  with  the  largest  Ecut.  We  hypothesize  the 
error  in  the  lattice  vectors  is  due  to  a  lack  of  van  der 
Waals  forces  for  non-overlapping  densities  in  current 
DFT  functionals.  This  is  supported  by  the  fact  that  the 
molecular  structures  are  generally  in  much  better 
agreement  with  experiment  than  are  the  unit  cell  vectors, 
which  should  depend  more  on  the  effect  of  dispersion 
forces  than  would  the  structure  of  an  isolated  molecule. 
Inaccurate  van  der  Waals  forces  cause  an  underestimation 
in  the  attractive  forces  between  the  molecules  in  the 
crystals  and  thus  an  overestimation  in  the  overall  lattice 
vectors.  This  deficiency  could  have  unforeseen 
consequences  on  all  molecular  crystal  calculations, 
particularly  for  those  at  low  pressure,  and  therefore 
caution  should  be  employed  whenever  interpreting  results 
obtained  from  any  solid  state  DFT  code.  To  properly 
describe  the  electronic  structure  of  the  crystals,  possible 
solutions  are:  (1)  new  functionals  need  to  be  developed 
which  include  accurate  dispersion  forces  for  non-  or 
nearly-overlapping  densities  (Rydberg  et  al.,  2003;  Dion 
et  al.,  2004),  (2)  employ  linear  response  theory  to 
calculate  the  frequency  dependent  susceptibility  and  use  it 
to  obtain  the  dispersion  forces,  (3)  or  use  perturbation 
theory  for  the  intermolecular  correlation  and  DFT  for  the 
intramolecular  correlation,  taking  care  to  avoid  double¬ 
counting. 
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